AD-A234  675 


NASA  Contractor  Report  187536 
1CASE  Report  No.  91-28 


ICASE 


ANALYSIS  AND  CONVERGENCE  OF  THE  MAC  SCHEME. 
1.  THE  LINEAR  PROBLEM 


R.  A.  Nicolaides 


Contract  No.  NAS  1-1 8605 
March  1991 


Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  Virginia  23665-5225 

Operated  by  the  Universities  Space  Research  Association 


NASA 

National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 

Hampton,  Virginia  23665  5225 


ANALYSIS  AND  CONVERGENCE  OF  THE  MAC  SCHEME. 

1.  THE  LINEAR  PROBLEM 


R.  A.  Nicolaides1 
Department  of  Mathematics 
Carnegie  Mellon  University 
Pittsburgh,  PA  15213 


ABSTRACT 

The  MAC  discretization  of  fluid  flow  is  analyzed  for  the  stationary  Stokes  equations.  It 
is  proved  that  the  discrete  approximations  do  in  fact  converge  to  the  exact  solutions  of  the 
flow  equations.  Estimates  using  mesh  dependent  norms  analogous  to  the  standard  H1  and 
L 2  norms  are  given  for  the  velocity  and  pressure  respectively. 
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1.  Introduction 


The  original  MAC  (Marker  and  Cell)  scheme  of  Harlow  and  Welch  is  now  over  thirty 
years  old.  Many  evolutionary  changes  have  occurred  and  the  MAC  approach  is  still  a  widely 
used  finite  difference  method  for  incompressible  flow  computation. 

The  MAC  scheme  is  described  in  [Hirt  1979]  and  in  most  CFD  books,  for  example 
[Fletcher  1989]  and  [Peyret  and  Taylor  1983].  It  also  forms  the  basis  of  several  flow  packages 
including  SOLA  [Hirt  1979]  and  SIMPLE  [Patankar  1980]. 

More  recently,  it  is  finite  element  techniques  which  have  dominated  incompressible  CFD. 
Although  highly  successful  in  many  ways,  the  absence  of  fully  satisfactory  lower  order 
schemes  continues  to  be  a  stimulus  to  the  development  of  new  approaches.  Recent  work 
on  finite  elements  may  be  found  in  [Girault  and  Raviart  1986],  [Gunzburger  1989]  and 
[Pironneau  1989]. 

The  stability  problems  of  low  order  finite  elements  [e.g  Boland  and  Nicolaides  1985] 
suggest  looking  for  stable  schemes  outside  the  finite  element  framework.  The  MAC  scheme 
is  one  possibility.  The  standard  MAC  method  suffers  from  being  limited  to  rectangular 
meshes.  A  second  difficulty  is  the  apparent  absence  of  any  rigorous  analysis  for  the  stationary 
problem  -  the  subject  of  this  report.  ([Porsching  1979]  discusses  the  evolutionary  problem.) 

On  the  first  point,  a  generalization  of  the  MAC  ideas  to  triangular  and  other  more 
general  meshes  was  presented  in  [Nicolaides  1989].  In  this  complementary  volume  or  covolume 
framework  the  MAC  ideas  are  given  a  control  volume  interpretation.  Complementary  pairs 
of  control  volumes  are  used,  and  they  can  conveniently  be  taken  as  the  triangles  and  polygons 
of  a  Delaunay- Voronoi  mesh  system.  When  specialized  to  a  standard  triangulation  of,  for 
example,  a  rectangular  domain  the  precise  MAC  equations  are  recovered.  [Nicolaides  and 
Wu  1991]  contains  an  implementation  for  a  specific  problem.  [Hall,  Cavendish,  and  Frey 
1991]  contains  a  more  general  implementation. 

The  covolume  formulation  is  based  on  discretization  of  the  vector  operators  grad,  div, 
and  curl  and  unlike  the  usual  MAC  derivation  is  coordinate  free.  The  error  estimates  that 
are  proved  here  follow  from  this  formulation  of  the  MAC  scheme.  A  related  approach  was 
used  in  [Nicolaides  1991]  to  discretize  and  estimate  some  div-curl  systems  and  several  results 
from  there  will  be  used  below. 

Before  proceeding,  a  few  technical  comments  are  in  order.  First,  we  will  present  the 
analysis  for  two  dimensional  problems  alt' lough  the  methods  used  are  not  inherently  two 
dimensional.  Second,  the  estimates  for  the  velocity  are  in  a  discrete  H1  norm  on  Cl.  This 
implies  in  particular,  estimates  for  the  vorticity  on  the  boundary.  While  such  estimates  are 
desirable  from  a  practical  point  of  view,  we  have  assumed  that  u  £  H2  and  p  6  H2 .  This  is 
more  regularity  than  one  would  like  to  see.  In  our  analysis,  the  MAC  scheme  itself  is  obtained 
as  a  discretization  of  the  “strong”  second  order  form  of  the  equations.  This  already  suggests 
the  need  for  extra  regularity.  The  third  point  is  that  while  we  use  constant  mesh  spacings 
in  the  x  and  y  directions,  the  domain  need  not  be  rectangular  or  convex.  Coordinate  free 
methods  make  this  possible.  Summation  by  parts  formulas  for  example  become  inherently 
multidimensional  and  not  limited  to  summations  along  grid  lines. 
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2.  Preliminaries 

For  clarity,  we  assume  a  rectangular  domain  f l  and  a  primal  Cartesian  mesh  having  x 
and  y  spacings  equal  to  h  and  h'.  To  avoid  unnecessary  complications,  we  will  suppose 
that  h  —  h' .  It  is  easy  to  modify  the  results  to  cover  the  contrary  case.  In  addition,  by 
joining  the  centers  of  adjacent  mesh  rectangles  (cells)  a  dual  (“staggered”)  mesh  is  made. 
At  boundaries,  nodes  of  the  dual  mesh  are  joined  to  the  midpoints  of  the  adjacent  boundary 
edges. 

The  N  nodes  of  the  primal  mesh  are  numbered  from  left  to  right,  and  from  bottom  to 
top.  The  T  nodes  of  the  dual  mesh  are  numbered  in  a  similar  way.  There  are  T  mesh  cells  in 
all,  and  there  are  N  dual  mesh  cells  including  the  half  size  cells  at  the  boundary.  Similarly, 
the  E  edges  of  the  primal  mesh  are  labelled  in  some  convenient  way.  There  is  a  bijective 
correspondence  between  dual  and  primal  edges,  each  edge  crossing  exactly  one  edge  of  the 
other  set.  The  numbering  of  the  dual  edges  reflects  this,  each  being  numbered  the  same  as  its 
primal  companion.  The  cells,  edges  and  nodes  of  the  primal  mesh  are  denoted  generically  by 
T;,  g j  and  v>ic  respectively.  Those  of  the  dual  mesh  are  similarly  denoted  by  primed  quantities 
such  as  a'y  The  lengths  of  the  dual  mesh  edges  are  equal  to  h  except  near  boundaries  where 
they  may  be  h/2.  A  direction  is  assigned  to  each  primal  edge  according  to  the  rule  that 
positive  is  from  low  to  high  node  number.  The  dual  edges  are  directed  by  the  convention 
that  (oj.,  a*.)  are  oriented  like  the  (x,  y )  axes  of  the  coordinate  system. 

Defined  on  each  primal  edge  Gj  is  a  component  of  the  velocity  field,  directed  in  the  positive 
direction  of  <r' .  This  velocity  component  is  computed  as  n  •  u  or  its  average,  where  u  denotes 
the  velocity  vector  at  the  midpoint  of  the  primal  edge,  and  n  is  the  normal  to  the  edge. 
Components  of  other  vector  fields  are  defined  similarly.  Such  sets  of  normal  components 
defined  on  the  primal  edges  can  be  identified  with  RE .  We  will  introduce  an  inner  product 
into  RE  by 

{u,v)w  :=  £  uivjWj. 

In  this,  Wj  equals  twice  the  area  of  the  figure  obtained  by  joining  the  nodes  of  a  primal  edge 
to  the  adjacent  dual  nodes.  The  sum  is  to  be  taken  over  all  edges  of  the  primal  mesh.  The 
associated  norm  is  denoted  by  ||  •  ||w  Clearly,  it  is  twice  a  discrete  L 2  norm.  This  inner 
product  space  is  denoted  by  U,  or  by  UQ  if  the  normal  components  assigned  to  the  boundary 
edges  are  all  zero. 

Various  scalar  fields  including  pressures  are  defined  at  the  dual  nodes.  They  can  be 
identified  with  elements  of  RT .  An  inner  product  on  RT  is  defined  by 

(<t>,  Q)a  :=  ^2  WiAi- 

Here,  the  sum  is  over  all  mesh  cells  and  A,  denotes  the  area  of  the  ith  cell.  The  associated 
norm  is  denoted  by  ||  •  || >\ .  This  inner  product  space  is  denoted  by  V. 

Analogous  to  V,  scalar  fields  defined  on  the  primal  nodes  can  be  identified  with  elements 
of  Rn  ,  and  an  inner  product  defined  by 

(V>,  x)a>  :=  £  fa'Xi'Ki 

r'€fi 
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where  the  sum  is  over  all  dual  cells  including  dual  boundary  cells,  and  A'k  denotes  the  area 
of  the  kth  dual  cell.  The  norm  is  denoted  by  ||  ■  jl^/,  and  the  inner  product  space  by  S,  or 
by  S0  if  the  boundary  values  are  all  zero.  Sometimes,  it  will  be  convenient  to  extend  this 
norm  and  inner  product  over  the  interior  or  boundary  nodes  separately.  A  subscript  will  be 
used  to  denote  such  dependence  thus:  ||  •  iU\n- 

For  each  primal  cell  rt,  discrete  flux  and  divergence  operators  are  defined  on  IA  by 

(Du)t  :=  Y  uj~h] 

<T}edTi 

and 

( Du\  =  (Du),/^. 

By  hj  v«  mean  hj  negatively  signed  if  the  corresponding  velocity  component  is  directed 
towards  che  inside  of  rt,  and  positively  signed  otherwise. 

For  each  interior  dual  cell  rj  discrete  circulation  and  curl  operators  are  defined  by 

(Cu),  :=  Y  ukh'k 

*‘kedr; 


(Cu),  =  (duWA'j. 

This  time  the  tilde  produces  a  negative  sign  if  the  corresponding  dual  edge  is  directed  against 
the  positive  sense  of  description  of  and  a  positive  sign  otherwise. 

This  definition  cannot  be  used  in  boundary  dual  cells.  To  extend  it,  we  must  require 
that  tangential  velocity  components  are  specified  along  boundary  segments  defined  by  the 
intersections  of  consecutive  dual  mesh  edges  with  I\  Then  we  can  define  discrete  circulations 
and  curls  in  the  same  way  even  for  the  boundary  dual  cells.  These  extensions  of  C  and  C 
to  the  boundary  are  denoted  by  Cj,  and  Cb-  We  will  consider  that  the  components  of  Cu 
are  associated  with  interior  nodes  V{  and  that  the  components  of  Ct,u  are  associated  with 
interior  or  boundary  nodes  as  appropriate. 

Normal  boundary  components  of  u  G  U  are  denoted  by  zdr  Sometimes,  it  will  be 
convenient  to  use  the  same  notation  to  include  tangential  components  ir  ne  sense  of  the 
previous  paragraph  as  well  as  normal  components.  It  will  be  explicitly  st  .ted  whenever  such 
tangential  boundary  components  are  included. 

Slope  operators  are  defined  on  V  and  S  by 


<f>l  J  -  <t>x 


where  ij  and  i2  are  the  nodes  defining  a[  and  the  positive  direction  is  from  ix  to  i2,  and 

W).  :=  11,  „  >  „ 

where  at  is  a  primal  edge  with  endpoints  and  z2. 

We  will  use  the  summation  by  parts  formulas  which  follow: 

(z)  (u,  Rip)w  =  (Cfctt,  Vu  G  U  Vip  e  S 


;i 


where  tangential  components  of  u  are  zero.  This  notation  is  consistent  with  the  definition  of 
Cj,  only  if  the  tangential  components  of  u  are  zero.  We  will  require  a  variant  of  this  to  cover 
the  case  when  Rip  is  restricted  to  interior  edges  -  those  with  at  most  one  node  on  T.  This  is 

(i)'  (u,  Rip)w  =  ( Cbu ,  ip)A,  Vu  G  Uq  Vt/>  €  S 

so  that  the  normal  boundary  values  of  u  are  zero  as  well. 

(u)  («,  G<p)w  =  -(Du,  <p)A  Vu  e  u0  vv  e  v. 

These  are  easily  proved  by  direct  computation. 


3.  Div-curl  Results 

A  covolume  method  and  analysis  was  given  in  [Nicolaides  1991]  for  the  div-curl  system 

divu  =  p 
curl  u  =  uj 

u  •  n|r  =  /. 

We  will  quote  some  results  of  [Nicolaides  1991]  which  are  used  below. 

In  the  notation  introduced  in  section  2  the  discrete  div-curl  approximation  is 

Du  ~  p 
Cu  —  U) 
u|r  =  / 

where  in  each  case  the  data  is  computed  by  simple  averaging  over  the  appropriate  geometrical 
element.  For  example, 

Q,  :=  — -  J  u  dx  dy 


fj  ■=  rjr  /  fds 

|  <7j  |  Jo  j 


where  <jj  is  numbered  anticlockwise  around  T. 

It  is  worth  remembering  that  Cu  are  the  discrete  curls  (normalized  circulations)  taken 
around  interior  nodes  only. 

If  fi  is  multiply  connected  there  is  an  additional  condition  [Nicolaides  1991],  but  we  are 
considering  rectangular  domains  here.  The  results  stated  below  are  valid  generally. 

Define 


:=  \k\  L u  ■ 


ids  a*  6  Q 


where  the  integration  is  along  the  positive  direction  t  of  a'k.  The  superscript  is  chosen  for 
consistency  with  [Nicolaides  1991].  If  ak  £  T  then  we  define 


u(2)  -  K 

uk  ■—  Jk- 
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An  error  estimate  follows: 


Theorem  1.  Assume  that  the  div-curl  equations  have  a  unique  solution  u  with  u  E  H2( Q). 
Then  (a)  the  discrete  div-curl  equations  have  unique  solution  u,  and  ( b )  the  estimate 

||u  -  u(2)||iv  <  Kh2 |u|Hj(n) 

holds  where  K  is  independent  of  u  and  h. 

Proof.  This  is  proved  in  the  remarks  following  Theorem  6.1  in  [Nicolaides  1991]. 


4.  Stokes  Equations 

The  inhomogeneous  stationary  Stokes  equations  are 

Au  —  yp  =  f  in  Q 

div  u  =0  infi 

ulr  =  g. 

Here,  17  denotes  a  polygonal  domain  in  R2,  and  T  denotes  its  boundary. 

The  standard  weak  form  of  the  Stokes  equations  has  unique  solution  u  €  H1(D),  p  E 
L2(0)  if  f  €  L2(Q)  and  g  E  H1,/2(r)  [Girault  and  Raviart  1986]. 

The  vorticity  u>  is  defined  to  be  curl  u  :=  dxv  —  dvu  where  the  velocity  field  u  :=  K  v)- 
In  terms  of  u>  and  using  the  incompressibility  constraint,  the  momentum  equation  becomes 

-curl  ui  -  VP  =  f  (2) 
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where  curl  denotes  the  operator  ( dv ,  —dx). 

Let  a,  denote  an  interior  mesh  segment,  and  let  t  and  n  be  the  unit  tangent  and  unit 
normal  directed  along  the  positive  directions  of  at  and  a\  respectively,  a i  may  have  at  most 
one  node  on  T.  Let  A  and  B  denote  ct.’s  nodes  where  the  positive  direction  is  from  A  to  B. 
Tal.mg  the  inner  product  of  the  last  equation  with  n  and  integrating  along  a,  from  A  to  B 
gives 

-w(J3)  ta;(4)  =  jf  ^dt  +  jf  f  nds.  (3) 

If.  say,  B  is  a  boundary  node  we  will  use  this  formula  to  extend  the  vorticity  to  B.  The 
extended  value  is  well  defined  and  independent  of  the  particular  line  segment  through  B 
(or  the  point  A  on  it)  under  the  assumptions  u ;  £  H2{ Q),  SJp  £  f  £  H1(f2).  To 

prove  this,  denote  the  extension  to  B  just  defined  by  uq(I?),  and  let  another  one  based  on  a 
positively  directed  line  segment  from  C  to  B  be  called  <jJi(B).  Then  it  follows  that 

Wi)-^(J)]  +  M>)-fl=/  (l^  +  f  n )dt-  I  (fUf-n )dt 

Jab  on  Jcb  on 

so  that 

M-B)  -  =  X(^flC)(Vp_f)'  nds 

=  0 

where  now  n  denotes  the  outer  unit  normal  to  the  boundary  d(ABC)  of  the  triangle  ABC 
and  the  trace  theorem  was  used  in  association  with  (2). 

Note  that  since  oj  £  by  Sobolev’s  lemma  u>  is  uniformly  continuous  on  0,  and  so 

can  be  uniquely  extended  to  be  continuous  on  fi.  The  extension  of  uj  onto  Q  using  (3)  and 
the  Sobolev  extension  coincide:  for  denoting  by  u ’s(B)  the  Sobolev  value  it  follows  that  with 
A  £  ft 

\v(B)  -  urs(B)\  <  |u>(  3)  -  o>(.4)|  +  \u>s{B)  -  w(.4)| 

<  f  \jr-\ds  +  |ws(B)  -  w.s(>4)| 

J  AS  UTl 

<  |A5|1/2!|p||H3(n)  +  ks(£)  -  ws(^)l 

where  the  integration  is  performed  along  the  line  segment  joining  A  and  B,  and  Cauchy’s 
inequality  and  a  trace  theorem  were  used.  The  result  follows  from  this. 

5.  Discrete  Stokes  Equations 

The  discretization  of  the  Stokes  equations  is  based  on  (3)  obtained  above.  The  momentum 
equation  may  be  written  as 

where 
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and  n  denotes  the  unit  normal  to  <r, ,  and  F  denotes  the  vector  of  average  normal  components 
of  f.  In  terms  of  the  operators  of  section  2,  the  discrete  equations  are  by  definition 

-Ru'-Gp'  =  F  a,  <=  fi 
Du1  —  0 
Cbu'  =  U)' . 

Boundary  values  for  u  are  defined  as  simple  averages  of  the  normal  and  tangential  com¬ 
ponents  of  g  along  the  boundary  edges  of  the  mesh,  so  that  a  typical  normal  boundary  value 
is 

u1  =  -  /  g  •  n  ds  a,  £  T 
h  J a; 

where  n  denotes  the  unit  normal  directed  positively  for  the  edge  and  h  denotes  the  edge 
length.  The  typical  tangential  boundary  value,  associated  to  a  primal  boundary  node  i\  is 
similarly  the  mean  of  the  tangential  component  of  g  between  the  dual  edge  intersections 
enclosing  i/j. 

To  normalize  the  pressure  approximation,  we  may  impose 


=  0. 

The  discrete  momentum  equation  has  one  component  per  interior  edge,  and  there  is 
one  incompressibility  constraint  for  each  primal  cell.  Since  there  is  one  unknown  velocity 
component  per  interior  edge  and  one  pressure  variable  per  cell  there  is  a  match  between 
equations  and  unknowns,  excluding  the  pressure  normalization. 

It  can  be  proved  by  direct  calculation  [Nicolaides  1988]  that  this  approximation  to  the 
Stokes  operator  is  precisely  the  MAC  approximation.  The  equations  would  be  identical  were 
it  not  for  our  treatment  of  the  data  by  exact  integration  in  place  of  the  quadrature  likely  to 
be  used  in  reality. 

Explicit  in  the  equations  above  is  an  unexpected  connection  with  the  velocity-vorticity 
formulation  of  the  Stokes  equations.  From  our  discrete  equations  above  we  see  that  the 
MAC  scheme  can  reasonably  be  described  as  a  discretization  of  that  formulation.  But  in  the 
standard  finite  difference  way,  they  are  derived  as  an  approximation  to  the  primitive  variable 
equations.  So  the  MAC  scheme  simultaneously  discretizes  the  two  forms  of  the  governing 
equations.  Whatever  advantages  or  disadvantages  are  perceived  in  either  formulation  are 
therefore  present  in  the  discrete  scheme.  From  a  slightly  different  viewpoint  it  follows  that 
there  are  discrete  analogs  in  the  MAC  framework  (and  in  the  covolume  framework  in  general) 
of  the  transformations  which  enable  us  to  go  from  one  formulation  to  the  other  [Choudhury 
and  Nicolaides  1991]. 

Including  the  pressure  normalization  there  is  one  more  equation  than  unknowns.  This  is 
not  always  convenient.  We  can  avoid  it  by  subtracting  the  mean  pressure  in  the  momentum 
equation.  For  this,  let  e  G  RT  be  the  vector  with  ones  in  all  its  positions,  and  let  A  €  RT 
have  Aj  :=  Aj/\Sl\.  In  place  of  Gp  we  write  G(I  —  eAl)p.  The  discrete  Stokes  equations 

become 


~RChu  G{1  -  e.4')p'  -  F 
Du  --  0 


( 


together  with  the  boundary  equations.  Now  the  pressure  normalization  is  built  in  and  the 
number  of  equations  and  unknowns  and  equations  is  the  same. 

Theorem  1.  The  equations 


—RCbu'  —  G(I  —  eA*)?'  =  F 
Du'  =  0 

with  prescribed  normal  and  tangential  boundary  values  u|r  have  a  unique  solution. 

Proof.  Consider  the  homogeneous  equations 

-RCbu'  -G(I  -eiV  =  0 
Du'  =  0 
u'|r  =  0. 

Included  in  the  homogeneous  boundary  condition  are  both  the  normal  and  tangential  values. 
Taking  the  inner  product  with  u'  €  Uq  we  obtain 

( u RCbu')w  +  (u1,  G(I  —  eA^p'^w  —  0. 

Using  the  summation  formulas  (i)'  and  (ii)  reduces  this  to 

\\M\\9a'  =  o 

and  in  particular  Cu'  =  0.  By  Theorem  3.1  (a)  it  follows  that  u'  =  0.  Then  G{I  —  eAt)pt  —  0 
and  (/  —  eA^p'  is  constant  since  the  mesh  is  connected.  But  p'  can  differ  from  its  mean 
by  a  constant  only  if  the  constant  is  zero,  and  so  p'  =  0.  Uniqueness  follows,  and  since  the 
coefficient  matrix  is  square  so  does  existence. 


6.  Error  Bounds 

We  will  define 


(w„). 

1 

J  uj  dxdy 

Vi  G  fl. 

(4) 

By  Stokes’  theorem  we  also  have 

("*)i 

1 

"4- 

/  u  •  t  ds 
hr' 

v,  €  JI. 

(5) 

Associated  with  u>v  is  a  unique  discrete  velocity  field  v  defined  by 


Cv  =  U)v 
Dv  —  0 
u|p  =  u' 
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where  we  recall,  that  u'jr  is  defined  as  the  simple  average  of  the  normal  component  g  •  n  on 
each  boundary  edge. 

The  field  v  defined  in  this  way  has  no  tangential  components  on  T  associated  with  it. 
We  will  now  define  them  to  be  equal  to  those  of  u' .  It  is  convenient  to  continue  to  call  the 
augmented  field  v.  Recall  that  the  tangential  components  of  u'  on  T  are  defined  as  averages 
of  the  prescribed  tangential  velocity.  It  follows  that  the  tangential  components  of  u'  —  v  are 
all  zero. 

The  last  equations  are  the  standard  covolume  approximation  to  the  system 

curl  u  =  u 
div  u  =  0 
u  •  njr  =  g  n. 

The  difference  u'  —  v  satisfies  the  discrete  div-curl  equations 

D(u'  —  v)  =  0 

C{u'  —  v)  —  oj'  —  u> 

(u1  -u)|r  =  0 

and  by  Theorem  3.2,  it  follows  that 

||u;  —  u||w  <  K\\u)'  —  U>t)||i4/.  (6) 

In  addition  to  the  interior  values,  boundary  values  for  u>„  are  defined  by 

(&v)i  :=  (Cbv)i  Vi  G  T. 

Subtracting  the  exact  and  approximate  momentum  equations  and  introducing  p  G  V 
equal  to  the  mean  pressure  in  the  primal  mesh  cells  we  obtain 

-  R(u'  -  Cbv)  -  G(p'  -p)  =  - R{u>  -  u>v)  -  -  GP )  a  e  ^  (7) 

The  meaning  of  the  notation  a  G  is  that  the  equation  has  one  component  for  each  edge  of 
the  orimal  mesh  which  is  not  on  the  boundary  T. 

It  follows  that 

-(u'  -  v,  R(u>'  -  u>„) )w  -  ( u '  ~  v,  G(p'  -  p))w  = 

I  <9 p 

(u'  -  V,  R{ u>  -  U>v))w  -  (u'  -  u,  /*(— )  -  Gp)w  a  G  n 

and  using  both  of  the  boundary  conditions  on  u'  —  v  and  the  summation  formulas  (i)'  and 
(zi)  that 

-{Cb(u'  -  v),  J  -  Gjv)a>  +  ( D(u '  -  v),  p  -  p)A  = 

~{Cb{u'  -v),u-  u>v)A>  - 


(«'  “  v>  -  Gp)w. 


Thus, 


dp 

("'  -  u/  -  "„)x'  =  ("  -  "  -  "u)x'  +  (u#  -  t\  /i(— )  -  Gp)w 

from  which  it  follows  by  Cauchy’s  inequality  and  bounding  ||u'  —  ul|iv  as  in  (6)  that 


||";  -  "«IU'  <  ||"  -  -i  K\\p(~)  —  Gp\\w-  (8) 

In  this  equation,  the  A'  norms  are  extended  over  all  nodes  including  those  on  I\ 

We  will  define  a  new  field  u*  6  U  as  follows:  on  interior  dual  edges  it  is  given  by  (1),  and 
on  dual  boundary  edges  by  a  similar  average.  It  also  has  tangential  bound  ry  values  equal 
to  those  of  u'.  By  definition  of  Cj,  it  follows  from  (4)  and  (5)  extended  to  T  that 


CbU*  =:  "  = 


Vi  £  fl. 


u>  gives  the  mean  vorticity  in  each  dual  cell,  including  those  on  T.  Notice  that  for  interior 
dual  cells  and  "  coincide. 

Using  (8)  we  now  have  with  Ke(p)  :  =  ||/i(|£)  —  Gp\\w 


\W  ~  ^lU'  — 
< 

< 

< 

< 


||"'  -  w„|U<  +  ||w„  -  q\\a. 

II"  -  "vlU'  +  ||"  -  "vlU'  +  Ke(p) 

||"  -  "|U'  +  2||"  -  a>v|U'  +  Ke(p) 

II"  —  "|U*  +  2||"  —  "v|U',r  +  Ke(p) 

II"  —  "IU',n  +  ||"  —  "lU'.r  +  2||"  -  "vlU'.r  +  Ke(p). 


(9) 


The  terms  bounding  the  error  can  be  estimated  by  approximation  theory  and  Theorem 

1. 


6.  Estimates 

To  begin  with  we  will  estimate  ||"  —  "H^n  in  (9)-  For  this,  consider  the  linear  functional 

/■1/2  ,1/2 

/j(")  :=  "(0,0)  -  /  /  d(  drj 

J- 1/2  J-l/2 

on  H2(Qi),  where  Qi  denotes  the  square  —1/2  <  {  <  1/2,  —1/2  <  tj  <  1/2.  It  is  clear  that 
the  functional  is  bounded  and  is  zero  on  linear  functions.  It  follows  that 

Ki(")l  < 

Changing  the  variables  in  the  integrals  to  x  :=  £h  and  y  rjh  gives  eventually 

|("  -  ")>|  <  X/i|"|j^(rj)  rj  €  £2. 

Squaring  and  summing  over  all  r'  £  fi  gives 

II"  -  "IU',o  <  -^^2|"U3(n)- 


10 


For  h  sufficiently  small  (<  |o»|Hi(n)|/l^|ffJ(n)|)  it  follows  that 

Ilw  -6>IU',n  < 

which  is  the  form  we  will  use  below. 

For  the  corresponding  boundary  term  we  consider  the  cell  — 1/2  <  £  <  1/2,  0  <  rj  <  1/2. 
Then  denoting  this  by  Q2,  we  have 

The  second  term  can  be  estimated  by  the  same  method  used  above  for  Qi.  For  the  first 
term, 

\  J*  Uj,(0,v)dv\  <  *K(0,  011(0.  J) 

<  K(\lu\Hi(q7)  +  |w|h»(<j,)) 

using  the  trace  theorem.  Changing  variables  to  x  :=  (h  and  y  :=  rjh  gives 

I  /  4  o/w(0,  y)dy\  <  if(|w|/ri(Qa>fc))  +  ^Mh’(q2,k))- 

JO 

Squaring,  multiplying  by  h2  and  summing  over  the  boundary  nodes  and  incorporating  the 
second  boundary  term  yields 

llw  —  —  #2(^4|w|ff*(n)  'r  h2\u\ Hi(n)) 

so  that  for  h  sufficiently  small  we  have  the  estimate 


||w  -  w|U',r  <  Kh\u}\H^Uy 

We  shall  amalgamate  the  interior  and  boundary  estimates  to  get 


||w  -  wll^n  <  Kh\u;\HHo) 

<  K"/i|u|Wj(n) 

for  h  sufficiently  small. 

Next  we  will  estimate  ||u>„  —  d>|j^,  r.  The  contribution  to  this  sum  associated  with  a  fixed 
boundary  node  i  has  the  form  |(u>„)i  —  (u>),|a  y.  Suppressing  the  dependence  on  i  we  have 

|wv-wj  =  IK-v)i~-(u*-v)2^  +  (u*  -  v)3A|/ y 

<  |(u*  -  u)i  -  (u*  -  v)2\/h  +  |(u*  -  u)3|/~. 

In  this,  the  subscripts  refer  to  the  two  dual  edges  of  length  ^  and  the  interior  dual  edge  of 
length  h  which  are  used  to  compute  the  discrete  curl  around  the  boundary  node.  There  is 
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no  contribution  from  the  boundary  itself  because  the  tangential  component  of  u*  —  v  is  zero 
there.  The  term  (u*  —  v)i,  for  example,  is  the  difference  between  the  average  normal  velocity 
component  measured  along  the  normal  from  the  boundary  to  the  dual  node  of  the  boundary 
cell  in  question  and  the  average  of  the  (prescribed)  normal  velocity  component  measured 
along  the  boundary  edge  of  the  same  boundary  cell.  ( u *  —  v)2  denotes  the  similar  difference 
of  averages  on  the  opposite  side  of  the  node  i.  The  union  of  these  two  boundary  cells  is  the 
cell  Q2)h  which  was  encountered  in  the  previous  estimate. 

To  estimate  the  first  term,  we  note  that  the  functional  (u*  —  v)i  —  (u*  —  v)2  is  bounded 
on  and  it  can  be  checked  that  it  vanishes  on  linear  vector  fields.  Using  the  familiar 

argument  and  a  scale  change  we  obtain  the  estimate 

IK  -  u)i  -  K  -  vh\  <  ^Mtf2(Qsh). 

Summing  the  corresponding  terms  over  T  gives 

E  l(«*  -  «).  -  K  -  »)»iJ  <  xwiu|jp(11). 

r 

For  the  remaining  part,  using  Theorem  3.1  we  obtain 

ElK-»)l!  =  *-’ElK-o)l*a 

r  r 

<  ^“2]CIK  -vW 

n 

<  K2h~2hA\u\n-i(n) 

<  K2h2\u\H*(ny 

In  this,  the  sum  on  the  left  is  over  all  dual  edges  parallel  to  T  appearing  in  the  definition  of 
the  boundary  circulations  (in  fact,  the  ‘dual  boundary’  of  fi)  and  the  right  hand  sum  in  the 
second  line  is  over  interior  edges  of  the  primal  mesh.  Assembling  the  last  two  estimates  we 
obtain 

ll^v  -  ^lU'.r  <  A'/i|u|wj(n)- 

Estimation  of  the  pressure  term  follows  similar  lines.  The  basic  functional  for  this  case 

h(p)  :=  /,/!  dr,  -  ( f'2  f  pdf  dr,  -  I'*  ['pdf  dp). 

J- 1/2  0(  J-1/2J-1  J -1/2  JO 

/3(  )  is  bounded  on  H2(Q3),  where  Q3  denotes  the  rectangle  with  corners  at  (  —  1,  ±1/2)  and 
(l,±l/2),  and  is  zero  for  linear  functions.  It  follows  that 

\h{p)\  <  •K'Ip!hj(Qs)- 


Introducing  the  scale  changes  1  :=  ( h ,  y  :=  Tjh  we  obtain 

,  1  /  V)  j..  Pfc2  ~  Pfc  1  ,  ^  *r|_i 

h  L  sr- dy  “ — r-1 5 
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where  the  integration  is  along  the  positive  direction  of  <7*  and  Qo(h)  denotes  the  scaled  region 
with  ok  separating  the  primal  mesh  squares  rkj  and  Tkl.  A  similar  result  holds  when  ok  is 
horizontal.  Squaring  and  summing,  it  follows  that 

iwf£)  -  Gpik  <  Kh\p\„,m.  (10) 

These  estimates  give  the  first  part  of  the  main  result: 

Theorem  1.  Under  the  regularity  assumptions  u  £  H2(fi),  u>  £  H2(Q)  and  p  £  H2(Q),  the 
estimate 

||w'  -  u)\\a<  <  Kh(\p\H*(n)  +  |u|W2(n)) 
holds  for  all  h  sufficiently  small. 

Recalling  that  ||u||^0)  =  ||div u|||,2(n)  +  l!cur^ ulli5(n)  ^  f°H°ws  that  the  norm  of  the 
error  is  indeed  a  discrete  (mesh  dependent)  Ho(fi)  norm. 

8.  Pressure  Error 

We  begin  this  section  by  recalling  the  following  standard  result: 

Lemma.  The  equation 

divv  =  /  £  L20(n) 

has  a  solution  v  €  Hj(fi)  satisfying 

llvllH‘(n)  <  K\\f\\L7(n). 

Proofs  may  be  found  in  [Girault  and  Raviart  1986],  and  in  [Temam  1984], 

We  will  apply  this  result  with 

/  ■=  p  —  p'  e  Ll{n) 

where  the  right  side  denotes  the  piecewise  constant  function  with  these  values  in  each  primal 
mesh  square.  Clearly 

Up  -  p'IU’(n)  =  IIp-p'IU 

so  that 

l|v||n.(n)  <  K\\p-p'\\/i.  (11) 

Next,  introduce  €  Ho  defined  on  each  edge  of  the  mesh  as  the  mean  of  the  normal 
component  of  v  on  the  associated  primal  edge,  i.e. 

:=  7  /  v  ’  n  d-3  crk  eCi. 

h  Jark 
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Use  of  the  divergence  theorem  shows  that 


Dv U)  =  p  —  p'.  (12) 

In  addition,  we  have 

l|v(,)lk  <  *I|v||h.(0)  <  K\\f  -  p'lU.  (13) 

Only  the  first  inequality  is  new.  It  is  a  consequence  of  the  fact  that  the  functional 

fl/3 

B\  :=  /  ui(0,  y)dy 

J- 1/2 

is  bounded  on  H1((54),  where  Qt  denotes  the  square  with  corners  (±1/2,  0)  and  (0,  ±1/2). 
Then  changing  the  scale  to  h  in  both  coordinate  directions  and  summing  over  all  the  dual 
edges  gives  the  result. 

Similar  to  i/1)  we  will  need  v*  which  is  defined  on  the  dual  edges,  as  the  mean  of  the 
tangential  component  along  the  dual  edge  in  question: 

v*k  :=  7~7T  [y-nds  a'k€fi. 

\ffk\  J*' 

n  points  along  cr'  in  this.  We  will  define  tangential  components  for  v*  on  the  boundary 
segments  delineated  by  dual  mesh  lines.  These  tangential  values  are  defined  to  be  zero. 
Now  we  have  ^  ^ 

(CbV*Y  —  I  v  ■  t  dt  ~  —  I  curl  v  dx  dy 
r'jJa rj  Tj  Jr'. 

from  which  it  follows  that 

<  ll«rlv|li,(0)  <  WSr.(n)  <  K\\P-P%-  (H) 

We  will  need  an  estimate  for  To  obtain  it,  we  first  note  the  estimate 

IK  -  v(1)||*y  <  #  (15) 

The  technique  for  proving  this  is  given  in  [Nicolaides  1991,  Section  6).  We  will  omit  the 
details  here.  Then 


UO(1)IU'  <  ||C,(»‘ -  +  ||C,v*|U. 

<  fiK  -  «(*)||^  +  ||C»V*|U 

<  m-p'\u 


by  (15),  (11)  and  (14)  . 

Now  taking  the  inner  product  of  v W  with  the  basic  error  equation  gives 
(v^,  R(u '  -  w))w  +  G(p  —  p'))w  =  (v^\  —  Gp)w 


(16) 
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and  using  the  summation  formula  (**)  gives 

Ilf  -  pTa  <  lltV’IU-IK  -  MU'.n  +  lMl)lkM§J)  -  opllw. 

Then  using  (16)  and  (13)  we  finally  obtain 

IIP  -  p'lU  <  K{\V  -  "IU',0  +  M§£)  -  Op\\w 

Using  Theorem  6.1  and  the  approximation  theory  estimates  we  now  have  the  second  part 
of  the  main  result: 

Theorem  1.  Under  the  regularity  assumptions  u  £  H2(ft),  u>  G  H2(Cl)  and  p  €  i/2(fl),  the 
estimate 

II P1  ~  P\\a  <  ^^(|p|HJ(n)  +  Mtfa(n)) 
holds  for  all  h  sufficiently  small. 

Thus  the  MAC  scheme  is  first  order  accurate  for  the  vorticity  and  pressure.  It  is  not 
known  yet  whether  the  velocity  is  of  higher  order  accuracy  than  the  vorticity.  Computations 
suggest  that  it  is  one  order  greater.  [Nicolaides  and  Wu  1991]. 
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